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Superflow of Bose-Einstein condensate in an optical lattice is represented by a Bloch wave, a plane 
wave with periodic modulation of the amplitude. We review the theoretical results on the interaction 
effects in the energy dispersion of the Bloch waves and in the linear stability of such waves. For 
sufficiently strong repulsion between the atoms, the lowest Bloch band develops a loop at the edge of 
the Brillouin zone, with the dramatic consequence of a finite probability of Landau-Zener tunneling 
even in the limit of a vanishing external force. Superfluidity can exist in the central region of the 
Brillouin zone in the presence of a repulsive interaction, beyond which Landau instability takes 
place where the system can lower its energy by making transition into states with smaller Bloch 
wavenumbers. In the outer part of the region of Landau instability, the Bloch waves are also 
dynamically unstable in the sense that a small initial deviation grows exponentially in time. In the 
inner region of Landau instability, a Bloch wave is dynamically stable in the absence of persistent 
external perturbations. Experimental implications of our findings will be discussed. 



PACS numbers: 03.75.Fi, 05.30.Jp, 67.40.Db, 03.65.-w 



I. INTRODUCTION 

As one of the most remarkable macroscopic quantum 
phenomena, superfluidity has attracted enormous atten- 
tion since its discovery [ll Q, The realization of Bose- 
Einstein condensation in dilute atomic gases Q has pro- 
vided physicists a new fertile ground for exploring many 
aspects of this fascinating phenomenon, including fric- 
tionless current Q and vortices In particular, there 
has been increasing interests and efforts to study super- 
fluidity and related phenomena in the system of a Bose- 
Einstein condensate (BECHn an optical lattice, such as 
Landau-Zener tunneling IljlJiJlliJosephson effect [HI, 
and dynamical instabilityjl^, H3, 13 ■ Recently, quan- 
tum phase transition between superfluidity and Mott- 
insulator^^ was observed in such a system. 

Superflow of BEC in an optical lattice is represented 
by a Bloch wave, which can be regarded as a plane wave 
modulated by the periodic potential. Bloch waves and 
Bloch bands are the basic language and concepts in pe- 
riodic systems. They are also essential to the study of 
superfluidity and many of its related phenomena in this 
periodic BEC system. 

One interesting phenomenon is the tunneling between 
Bloch bands at the edge of the Brillouin zone. Previ- 
ously, people studied experimentally this inter-band tun- 
neling by dragging cold atoms with accelerated optical 
lattices The cold atoms are very dilute and the in- 
teraction between them can be ignored entirely. It is 
then curious to know how the tunneling will change if the 
cold atoms are replaced with a BEC, where the atomic 
density is relatively large and the interaction between 
atoms can no longer be neglected. Our study shows that 



the interaction will increase the tunneling probability in 
general p. HSj. When the repulsive interaction is over a 
critical value, a loop appears in the lowest Bloch band, 
resulting in non-zero tunne ling probability even in the 
adiabatic limit [alia liii I2ill2il< This breakdown of adi- 
abaticity is the result of superfluidity and can be viewed 
as a hysteresis phenomenon |3. 12^ . 

The most dramatic manifestation of superfluidity is 
that the boson system can maintain its current speed in 
a very tight space, such as a narrow capillary. Landau 
proposed a criterion for superfluidity [2|, y|: if the cur- 
rent moves slower than sound, it experiences no viscosity; 
otherwise, the system suffers an instability and loses its 
speed. In a homogeneous BEC, the sound speed is pro- 
portional to the square root of the BEC density. It is not 
clear how superfluidity can be defined and studied for a 
BEC in an optical lattice. Our method is to examine the 
energy dispersion of Bloch waves. When the Bloch waves 
are energy minima of the system, they represent super- 
flows; when the Bloch waves are energy saddle points, 
they suffer Landau instability: external disturbances can 
render the system to emit phonons thus reduce its speed. 

The BEC Bloch waves can be dynamically unstable, 
that is, the system diverges away from the original Bloch 
state upon small instantaneous perturbation|l3(. This 
instability does not exist either for a BEC in a free 
space or for a free particle in a periodic lattice; it hap- 
pens only when there are both interaction and periodic 
lattice. The dynamical instability has been observed 
experimentally ^2 E3- an d further confirmed by other 
theoretical studies and experimental observations plll24L 
|25| . This dynamical instability is quite a general phe- 
nomenon and widely exists in other systems, such as a 
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BEC confined in a ring[26l |27]]. Its general implication 
and relation to superfluidity will be discussed in this re- 
view. 

We will discuss these phenomena in detail and review 
many interesting results for the quasi-one dimensional 
case, where the optical lattice is created by two counter- 
propagating off-resonance lasers while the lateral motion 
of the BEC can be either ignored or confined [7I HoL l2Sj . 
In Section II, a brief introduction of the mean-field the- 
ory of BEC systems is given for the sake of completeness 
and introduction of the notations. In Section III, we de- 
fine Bloch waves and Bloch bands for BECs in optical 
lattices and present some general results and the numer- 
ical methods to find them. In Section IV, we simplify 
the system to a two-level model to study the tunneling 
between Bloch bands. We point out the tunneling in the 
adiabatic limit is related to a general problem, adiabatic 
evolution of nonlinear quantum states. In Section V, the 
Landau instability and dynamical instability of the Bloch 
waves are studied. Other approaches to superfluidity and 
experimental observation and general implication of dy- 
namical instability are discussed. Finally in Section VI, 
we summarize the paper. 



II. MEAN-FIELD THEORY OF BEC SYSTEMS 

Even though the BEC is quite dense compared to the 
cold atoms, it is still very dilute with typical densities at 
10 -20 m -3 , which is thousands of times more dilute than 
the air. Due to this diluteness, along with the extreme 
low temperature, the interaction between atoms can be 
characterized by the s-wave scattering and modeled with 
a (5-function. When the number of atoms is large, the 
BEC system is very well described by the mean-field the- 
ory given by the celebrated Gross-Pitaevskii eauationp^ 

.^dtb h 2 d 2 ip s , 4:7rh 2 a s . ,. 0l ,„ . 

in ^ = -n=ir5 + v (*)* + -H 2 V>, (2.1) 

at 2m ox m 

where m is the atomic mass, a s is the s-wave scattering 
length, and V(x) is the external potential imposed on the 
system. In our case, the potential is the optical lattice 
created by two laser beams and is given by 



V(x) = V cos(2k L x) , 



(2.2) 



where Vq is proportional to the laser intensity and k^ is 
the wave number of the laser. 

For simplicity, we will instead use the following dimen- 
sionless Gross-Pitaevskii (CP) equation 



.dip 



1 d 2 ip 



(2.3) 



In the above equation, all the variables are scaled to be 
dimensionless with the system's basic parameters: the 

strength of the periodic potential v is in units of -, 

the wave function ip in units of Juq where no is the 



averaged BEC density, x in units of ^j— , t in units of 
j^r. The coupling constant c = ^g^. L 

This system can also be regarded as a Hamiltonian 
system governed by the grand canonical Hamiltonian 

H=Jdx{r(- \^+vcos{x))^+ c -w - m 2 }, 

where fj, is the chemical potential. The GP equation 
(|2.3() can be obtained by variation of the Hamiltonian, 
idip/dt = SH/Sip*. 



III. BLOCH WAVES AND BLOCH BANDS 

Among all possible solutions of Eci. (|2.3|) . there are 
states which extremize the Hamiltonian (|2.4[1 and hence 
satisfy the time-independent Gross-Pitaevskii equation 



ld 2 1p 1 ,,2 , 



v cos(x) ip = (lip . 



(3.1) 



We call these extremum states nonlinear eigen- 
states, which are also called nonlinear coherent modes 
elsewhere |30|: accordingly, \x are nonlinear eigenvalues. 
Bloch waves are the nonlinear eigenstates of the follow- 
ing form 



ip(x) 



ikx 



4>k(x) , 



(3.2) 



where 4> k (%) is a periodic function of period 2-7T and k is 
the Bloch wave number. In particular, from Eq. l|3.1|l we 
have the following equation for each Bloch wave state <j>k 

- \(-r + ik ) 2( Pk + c\(p k \ 2 (p k +vcos(x) cf) k = fJ-(k) 4> k ■ 
2 ax 

(3.3) 

The set of eigenvalues fi(k) then form Bloch bands. 

When v = 0, the plane waves are the solutions of 
Ea. (|3.1() . and they represent a BEC flow with speed of 
k. As is well known, when the speed is smaller than the 
speed of sound, k < y/c, it is a superflow; when k > yfc, it 
develops a Landau instability and loses the superfluidity. 
Bloch waves are the counterpart of these plane waves in 
a periodic system. Whether these Bloch waves represent 
superflows or not is determined by the energy dispersion 
of these Bloch waves as we will show later. 

In the linear case, c = 0, all eigenstates are Bloch 
waves [3lj. The situation is very different for the periodic 
nonlinear system, where it is possible to have non-Bloch 
wave eigenstates. Furthermore, there are nonlinear Bloch 
waves that do not have their linear counterparts. Here 
is a beautiful example |ljj. Eq.ETJ has an exact 
solution of Bloch wave at the Brillouin zone edge k = 1/2, 



1> 



\Jc + v + \fc- 



y/c + V — yfc- 



(3.4) 

This Bloch wave state only exists when c > v so that 
it has no linear counterpart. This shows that there are 
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some "extra" nonlinear Bloch waves when the nonlinear- 
ity is strong enough. Their corresponding "extra" eigen- 
values should make the nonlinear Bloch band look very 
different from the linear Bloch band. This is indeed the 
case: our numerical results show a loop structure in the 
lowest Bloch band when c > v as seen in FigQ This 
loop structure was first found in Ref . |a| with a two- mode 
approximation, then confirmed in Ref. 20|, |21| . In partic- 
ular, in Ref. |2lj . a loop is also found for the second band 
at the Brillouin zone center. 
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wave 4>k- The method described in Ref.[13j is to expand 
(fik in Fourier series and minimize the Hamiltonian 12. 4|) 
in the space expanded by Fourier coefficients. Another 
one is used in Ref. po| . where Ea. (|3.1|l is first solved for 
different values of /i and k then the Bloch band fi(k) is 
found by interpolation. 

A much better method is the following. One still starts 
with the expansion of <f>k in Fourier series, 

JV 

<t> k (x)= E a > e ""< ( 3 - 5 ) 

n=-JV 

where N is the cut-off. With the substitution of the above 
equation into Eq. (|3.3() , one obtains 2 N + 1 equalities for 
the coefficients of each Fourier term e mx , 

fn(ao,a±i, ■ ■ ■ ,a±N,fi) = 0. (3.6) 

Finally, the Bloch wave is found by minimizing the fol- 
lowing sum 

N 

S= £ f n - (3.7) 

n=-JV 

The results in Fig|T] are obtained with this method by 
using N — 10. 



FIG. 1: Lowest Bloch bands at v = 0.1 for c = 0.0, c = 0.05, 
c = 0.1, and c = 0.2 (from bottom to top). As c increases, 
the tip of the Bloch band turns from round to sharp at the 
critical value c = v, followed by the emergence of a loop. 

This loop structure is a consequence of superfluidity. 
One clear indication is that the exact Bloch wave at the 
Brillouin zone edge, Ea. (|3.4|l . represents a flow of non- 
zero speed, \f c 2 — v 2 /2c 19]. This is quite remarkable 
if we recall that all the Bloch waves at the zone edge 
carry no currents in the linear case. For a free particle, 
the flow e lx / 2 is completely stopped due to first-order 
Bragg scattering by the periodic potential, leading a zero- 
current Bloch waves at the edge, k = 1/2. However, when 
the superfluidity gets stronger, that is, when c > v, the 
flow e lx / 2 can no longer be stopped by Bragg scattering, 
yielding a zone-edge Bloch wave that carries a current. 
The connection between superfluidity and the loop struc- 
ture is discussed in-depth in the context of hysteresis by 
Mueller |3. 

These Bloch states can be prepared experimentally, at 
least for the ones in the lowest Bloch band by adiabatic 
control. For a trapped cigar-shaped BEC, we slowly turn 
on an optical lattice along its longitudinal direction; we 
then have a BEC in the Bloch state at the center of the 
Brillouin zone. Other Bloch states can be achieved by ac- 
celerating the lattice for a certain amount of time. These 
experimental techniques of adiabatic turning-on and ac- 
celerating optical lattices have been demonstrated suc- 
cessfully with either cold atoms [H HI or BECspjo). 

There are several numerical methods to find the Bloch 



IV. TUNNELING BETWEEN BLOCH BANDS 

Consider the scenario that we drag the BEC across 
the entire Brillouin zone with an accelerating lattice. If 
the BEC is initially in a Bloch state belonging to the 
lowest Bloch band and the acceleration is small enough, 
the system will stay in the band and undergo Bloch 
oscillations 

HEi- As one increases the acceleration, the 
BEC will have increasing chance tunneling into the upper 
band at the edge of the Brillouin zone, where the band 
gap is the smallest. Without interaction, this tunneling 
is nothing but the famous Landau-Zener tunneling 33] , 
a well understood phenomenon. With the interaction in 
our case, the tunneling behavior is strongly modified, in 
particular, by the loop structure in the Bloch band as we 
will see later. 

In an accelerating lattice, the system is described by 
d 1 d 

i^ip = --(-zr + iat) 2 ip + v cos(x)^ + c\ip\ 2 ip , (4.1) 
ot 2 ox 

where a is the acceleration. Even in the linear case, c = 0, 
the above equation is difficult to analyze. To reduce the 
mathematical complexity, we will approximate it with a 
two-level model without losing the essential physics. 

A. Two-level Model 

The tunneling mainly occurs around the edge of the 
Brillouin zone, k = 1/2, where the wave function is dom- 
inated by two plane wave components. We zoom in on 
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this region and write 

cj>(x, t) = a{t)e ikx + b{t)e l(k ~ 1)x , (4.2) 

where \a\ 2 + \b\ 2 = 1. Following the simple algebra de- 
scribed in Ref.jl] by substituting the above two- mode 
wave function into Ea. (|4.1() . we obtain a two-level non- 
linear Schrodinger equation 

'1(0-^(0 ■ (4 - 3) 

where 

"w-K 7+ . M -,-«)■ (4 - 4) 

where 7 = ai and k = |6| 2 — |o| 2 . 

However simplified it may appear, this two-level model 
captures the essence of our BEC system as it reproduces 
the looped Bloch band (see FigJ2J . In this simple model, 
the "Bloch bands" are obtained by finding the eigenval- 
ues of H(j), 

#(7)(fc) =M(7)(j) • (4-5) 

The eigenvalues ^(7) are also called adiabatic energy lev- 
els. How the levels change with the coupling strength c 
is shown in Fig|2 where we again see the appearance of 
the loop for c > v in the lower energy level. 
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FIG. 2: Adiabatic energy levels for c = 0.05 (left) and for 
c = 0.2 (right) at v = 0.1. For comparison, the levels for the 
linear case, c = 0.0, are plotted in dashed lines. 

For this two-level model, the original tunneling prob- 
lem in the BEC system can be restated as follows: if 
the system starts in the lower level, as the energy bias 7 
changes with the sweeping rate a from the far negative 
end to the far positive end, what is probability that the 
system ends up in the upper level? 



When c = 0, this is precisely the well-known Landau- 
Zener model, which can be solved exactly. The tunneling 
probability is[33| 

r =exp-— -. 4.6 
2a 

It is clear from this formula that, in the adiabatic limit 
a — * 0, the tunneling probability is zero, that is, the sys- 
tem stays in the lower level throughout the entire pro- 
cess. This is nothing but a special case of the quantum 
adiabatic theorem|34|: in the adiabatic change of a pa- 
rameter the system starting in an eigenstate stays in the 
same instantaneous eigenstate. 

When c > v, this adiabaticity is broken down by the 
loop structure: Suppose we start with a state on the 
lower adiabatic level, and move it up along the branch 
by changing 7 so slowly such that little tunneling to the 
upper level is generated. After passing the crossing point 
X, the state remains in the course moving up in energy 
until hitting the terminal point T, where it has no way 
to go any further except to jump to the upper and lower 
levels. This leads to a nonzero probability tunneling into 
the upper level in the adiabatic limit. 

This qualitative analysis is corroborated by our numer- 
ical integration of Ea. H4.3|) . The tunneling probabilities 
as a function of the sweeping rate a are plotted in Fig|3| 
for different values of c at v — 0.1. In general, the tunnel- 
ing probability is enhanced due to the interaction. How- 
ever, its dependence on a changes dramatically with the 
increasing interaction strength c. When c < v, the tun- 
neling probability approaches zero in the adiabatic limit, 
a — > 0. It appears that the approaching is exponential, as 
supported by a recent analysis [Tsj. At the critical point, 
c — v, the tunneling probability still turns to zero as the 
sweeping slows down, however, not exponentially. 

When c is over the critical value, it becomes apparent 
in Fig|3 that the tunneling curve tends to intersect the 
vertical axis at a finite value, indicating the tunneling 
probability is not zero at the adiabatic limit, a — * 0. Of 
course, one can never be absolutely sure that the tun- 
neling probability at a — > is not zero by numerical 
calculation since there is always a limit on the smallness 
of a in a computer code. This uncertainty is put to an 
end in Ref.|l8^. where this nonzero adiabatic tunneling 
probability is obtained analytically. 

One can also directly integrate numerically the orig- 
inal Schrodinger equation Ea. (|4.1|l to find out the tun- 
neling probability and how it depends on the accelera- 
tion a. This is done in Ref . |35| . where one sees a similar 
crossover from the exponential to non-exponential depen- 
dence of the tunneling probability on a. We argue that 
this crossover behavior can serve as an indirect evidence 
of the existence of the loop structure in the Bloch band. 
We believe that this crossover should be observable in an 
experiment similar to the one conducted in Ref.|l0j]. 

The nonlinear Landau-Zener tunneling studied here is 
likely as general as the linear counterpart and should 
exist in many physical systems. We found one exam- 
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FIG. 3: Numerical results of tunneling probability as a func- 
tion of q at v = 0.1 for different values of c. For comparison, 
the exact result Eg. 14.61 for the linear case c = is plotted 
in solid line. The small wiggles are not physical and they are 
resulted from the numerical errors of our computer codes. 



pie in molecular magnets |36(, where the jamming effect 
among interacting molecular spins is the result of nonlin- 
ear Landau-Zener tunneling. 



B. Nonlinear quantum adiabatic theorem 

The nonlinear two-level model l|4. 3JI is a very good ex- 
ample of a general problem: how a nonlinear quantum 
system, such as a BEC system, responses to the change 
of a system parameter. Of particular interest is that the 
parameter change is adiabatic since the adiabatic con- 
trol is an important tool in preparing a quantum system 
(linear or nonlinear) into a desired quantum state |37j . 
For a linear quantum system, the answer is the quantum 
adiabatic theorem, which can be found in a standard 
textbook j3^|. The theorem states that if the system is 
initially in a non-degenerate eigenstate, it stays in the 
same eigenstate as the parameter changes adiabatically. 

On the other hand, the answer for nonlinear quan- 
tum systems is not as straightforward. The nonlinear- 
ity proves to be a big obstacle for generalizing the linear 
quantum adiabatic theorem 38]. Very non-trivial mod- 
ifications to the linear case are expected as it is found 
in the nonlinear two-level model (|4.3|) that tunneling be- 
tween energy levels hap pens even in the adiabatic limit. 

In a recent paper|39|, we have successfully generalized 
the quantum adiabatic theorem to nonlinear quantum 
systems by combining ideas from classical adiabatic dy- 
namics and quantum geometric phases. It is found that 
the eigenstates correspond to fixed points of an equiva- 



lent classical problem. The adiabatic evolution of these 
eigenstates depends on the nature of the corresponding 
fixed points. If the fixed point is elliptic, the system can 
stay in the same eigenstate as in the linear case. Oth- 
erwise, when the fixed point is hyperbolic, the system 
will stray away from the eigenstate, leading to tunneling 
between different eigenstates. This adiabatic condition 
can also be specified in terms of Bogoliubov spectrum. If 
the Bogoliubov spectrum is real, the adiabaticity is kept; 
otherwise, it is broken. 



V. SUPERFLUIDITY AND DYNAMICAL 
INSTABILITY 

A. Superfluidity and Landau instability 

A quantum Bose liquid or gas at very low temper- 
atures possesses a very remarkable property known as 
superfluidity: its flows suffers no viscosity through cap- 
illaries or other types of tight spaces when its speed is 
slower than a critical value. This was discovered even 
before the conception of Bose-Einstein condensation p|. 
Landaud gave a simple explanation for this intrigu- 
ing phenomenon. He argued that, a quantum current 
suffers friction only when the system's elementary exci- 
tation, phonon, can lower its energy. With the Galilean 
transformation, he found that a quantum liquid flowing 
with a speed smaller than the sound speed enjoys such 
frictionless transport. When it flows with a speed larger 
than the sound speed, the quantum liquid can lower its 
energy by emitting phonons, leading to the slowdown of 
its speed. 

Such lose of superfluidity can also be regarded as an 
instability: the system can no longer stay in its original 
state thus keep its flow against perturbation of the rough- 
ness of capillaries or other external disturbance. We will 
call this instability Landau instability. Note that this 
Landau instability is different from the Landau instabil- 
ity discussed in plasma physics [4fj|. 

In our case, the study of superfluidity or Landau in- 
stability is complicated by the presence of the optical 
lattice, which breaks the translational symmetry along 
its flow direction and modulates the BEC density. This 
complication is only mathematical. The physics is still 
the same, that is, the Landau instability is determined 
by whether the elementary excitation around a BEC flow 
lowers the system's energy or not. If it does not, the BEC 
flow is a local energy minimum of the system and it is a 
superflow. Otherwise, the flow is an energy saddle point 
and it suffers Landau instability (see Fig0J . 

Whether these Bloch waves are local energy minima is 
determined by their energy dispersions. To calculate the 
dispersion, we perturb the system around a Bloch state 



ij>(x) = e ikx <t>h{x) + 8<t> k {x) 



(5.1) 



Due to the periodicity of our system and the Bloch 
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Energy local minimum Energy saddle point 



Dynamical Instability 




time time 



FIG. 4: Schematic illustration of superfluidity, Landau in- 
stability, and dynamical instability. A dynamically unstable 
state must be a saddle point; a state that is an energy saddle 
point can be dynamically stable. 

waves, these perturbation can be decomposed into dif- 
ferent modes labeled by q, 

8<t> k {x,q) =u k (x,q)e i '> x +v* k (x,q)e- i « x , (5.2) 

where q is also a kind of Bloch wave number and ranges 
between —1/2 and 1/2. The perturbation functions u k 
and v k are of periodicity of 2-7T in x. By substituting 
the perturbed state l|5.1|) into the Hamiltonian (|2.4|) and 
neglecting terms of orders higher than two, we obtain a 
quadratic form of the energy deviation, which is block 
diagonal in q. Each block is given by 

SE k = dx ( u* k , v* k ) M k (q) (fy, (5.3) 

where 

A> <5 ' 4) 

with 

1 d 

C(k) =--(— + ikf + v cos(x) -n+ 2c|0 fe | 2 . (5.5) 

2 ox 

If M k {q) is positive definite for all —1/2 < q < 1/2, 
the Bloch wave 4> k is a local minimum and represents 
a superflow . Otherwise, SE k can be negative for some 
q, indicating that the Bloch wave is a saddle point and 
describes only a normal flow. 

The special case v = 0, when the optical lattice is 
turned off, is the simple case considered by Landau[2j. 
In this case, the Bloch state <f> k becomes a plane wave 
e lkx and the operator M k (q) becomes a 2 x 2 matrix 

Mi(9) = ( V/2+* 5 +c ?2/2 _ Vc ) (5 . 6) 



whose eigenvalues are 



A ± = | + c±vW+^. (5.7) 

Since A+ is always positive, the matrix M k {q) is not 
positive definite only when A_ < 0, or equivalently, 
\k\> \/q 2 /4: + c. It immediately follows that the BEC 
flow e tkx becomes a saddle point when the flow speed ex- 
ceeds the sound speed, \k\ > y/c. This recovers exactly 
the Landau condition for the breakdown of superfluidity 
, which has recently been confirmed experimentally on 
BEC @. 

The situation becomes more complicated in the case of 
v > 0, when the optical lattice is turned on. We resort to 
numerical calculations to examine whether the matrices 
M k (q) are positive definite and the corresponding BEC 
Bloch waves are local energy minima. Our focus is on 
the Bloch states in the lowest Bloch band, excluding the 
loop. The results are summarized in the stability phase 
diagrams in FigEI where a wide range of values of v and 
c are considered in order to give a global picture of how 
the instabilities change with the two system parameters, 
c and v. The results have reflection symmetry in k and 
q, so we only show the parameter region, < k < 1/2 
and < q < 1/2. 

In the shaded area (light or dark) of each panel of Fig [3] 
the matrix M k (q) has negative eigenvalues, and the corre- 
sponding Bloch states <fi k are saddle points. For those val- 
ues of k outside the shaded area, the Bloch states are local 
energy minima and represent superflows. The superflow 
region expands with increasing atomic interaction c, and 
occupies the entire Brillouin zone for sufficiently large c. 
On the other hand, the lattice potential strength v does 
not affect the superflow region very much as we see in 
each row. The phase boundaries for v <C 1 are well re- 
produced from the analytical expression k = \J q 2 /4 + c 
for v = 0, which is plotted as triangles in the first column. 

The Landau instability of the states in the lower branch 
of the loop has been discussed in Ref . [2l| . It is found that 
the area of superfluidity increases when the interaction c 
increases or the lattice strength v decreases. 

B. Remarks on superfluidity 

So far, we have been identifying superfluidity as stabil- 
ity in the sense of Landau0, yj. In this approach, there 
is an implicit assumption that the boson system has al- 
ready achieved phase coherence and can be described 
by a complex order parameter, the macroscopic wave 
function. Whether the state described by the macro- 
scopic wave function possesses superfluidity is then ex- 
amined by how the system responses to perturbations. 
With the assumption of phase coherence, this approach is 
not well-equipped to study quantum phase transitions at 
zero temperature, such as the superfluid-Mott-insulator 
transition[l5l Il6|. Nevertheless, this approach can still 
provide some signals for the onset of some possible phase 
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v=0.01 



v=0.05 



v=0.10 




FIG. 5: Stability phase diagram of BEC Bloch waves, k is the 
wave number of BEC Bloch waves; q denotes the wave number 
of perturbation modes. In the shaded (light or dark) area, 
the perturbation mode has negative excitation energy; in the 
dark shaded area, the mode grows or decays exponentially 
in time. The triangles in (a.l-a.4) represent the boundary, 
q 2 /4 + c = k 2 , of saddle point regions at v = 0. The solid 
dots in the first column are from the analytical results of 
Eq. lEnSl . The circles in (b.l) and (c.l) are based on the 
analytical expression 15.141 . The dashed lines indicate the 
most unstable modes for each Bloch wave <j>k- 



transition as we will discuss in the context of dynamical 
instability in the next section. 

There are various other approaches to superfluidity. 
One of particular interest was proposed by Fisher et al. 
in Ref.|41|. where superfluidity is defined as how strong 
the system is affected by an artificially imposed phase- 
twist at boundaries. Since only the superfluid part of the 
system responses to the phase-twist, by calculating the 
change of the system energy, one can find out the super- 
fluid fraction of the Bose system. This approach has its 
root in the Josephson effect, where the phase coherence 
aspect of superfluidity was exposed for the first time[4^|. 

This phase-twist approach does not assume the phase 
coherence; to the opposite, it is particularly designed 
to calculate how much phase coherence a system has. 
It has been widely used in the Bose-Hubbard model to 



study the superfluid-Mott-insulator transition ^Sj. With 
the recognition that a BEC in an optical lattice can 
be described by the Bose-Hubbard model[l(j> this ap- 
proach has now been adopted to study the superfluid- 
Mott-insulator transition for various systems of BECs in 
optical lattices |4ll4o|. 

Both approaches are of perturbative nature, that is, 
they study how the Bose system responses to small per- 
turbations. In the Landau's approach, all possible dis- 
turbances are considered; one is able to determine the 
robustness of a superfluid as measured by the critical ve- 
locity. In contrast, a special type of perturbation is con- 
sidered in the method of phase- twist, namely, an artificial 
phase twist imposed at the boundaries. This perturba- 
tion is sensitive to the interplay between kinetic energy 
and interaction. It is not clear how this phase twisting 
corresponds to usual perturbations (such as disorder or 
roughness) which would degrade the flow of a normal 
fluid. In short, the relation between these two pictures is 
still a subject of active debate [46|. 



C. Dynamical instability 

One unique feature in the system of BECs in optical 
lattices is dynamical instability, which does not exist in 
the absence of either atomic interaction, c = 0, or an op- 
tical lattice, v = 0. Many of the BEC Bloch waves can be 
dynamically unstable against certain perturbation modes 
q only when both factors are present. By dynamical in- 
stability, we mean that small deviations from a state grow 
exponentially in the course of time evolution (see Fig0J. 

The dynamical instability can be determined from the 
linear stability analysis of the GP equation Eq. i|2.3|) . As- 
sume that the system experience a small disturbance 5(f>k 
at a Bloch state, 



ip{x,t) 



4> k (x) + 54> k (x,t) 



(5.8) 



where the disturbance can be similarly written as 

5(j> k (x,q)= u k (x, q, t)e^ x + v* k (x, q, tje"* 9 * . (5.9) 

Plugging the above wave function into Eq. I|2.3[) and keep- 
ing only the linear terms, we arrive at a linear dynamical 
equation, 



d 



where 



<JzM k {q) 



I 

-/ 



Uk 



(5.10) 



(5.11) 



The eigenvalues of the matrix o z M k {q) determine 

the dynamical instability: If they are real for all —1/2 < 
q < 1/2, the Bloch state is dynamically stable; if other- 
wise, that is, some of them are complex, it is dynamically 
unstable. 
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One who is familiar with the quasi-particle excitation 
in quantum fluid may immediately notice that a z Mk(q) 
also appears in the Bogoliubov approach and gives us 
the spectrum of the phonon excitation 2]. When all of 
its eigenvalues £k(q) are real, the positive half are the 
phonon spectrum while the other half are deemed non- 
physical. However, we will refer to the modes correspond- 
ing to the negative half as anti-phonons for the sake of 
easy reference in the following discussion. 

Again, we first look at the case v = 0, where the eigen- 
values of cr z M k (q) are 

£±(q) = kq± y/q2 C + q*/4. (5.12) 

These eigenvalues are always real, implying that the BEC 
flows in a free space are always dynamically stable. When 
v^O, the situation is totally different: the eigenvalues 
Ek{q) of a z Mk(q) can be complex and Bloch states can 
be dynamically unstable. The results are summarized in 
Fig.|5l where these Sk (q) are complex in the dark-shaded 
areas. The dark shaded areas lie completely inside the 
light shaded areas, as the result of a rigorous conclusion 
that only saddle-point Bloch states can have dynamical 
instability. Its proof is given in the appendix. 

The dynamical instability is the result of the resonance 
coupling between a phonon mode and an anti-phonon 
mode by first-order Bragg scattering. As shown in the 
appendix, the matrix a z Mk(q) is real in the momentum 
representation, meaning its complex eigenvalue can ap- 
pear only in conjugate pairs and they must come from a 
pair of real eigenvalues that are degenerate prior to the 
coupling. Degeneracies or resonances within the phonon 
spectrum or within the anti-phonon spectrum do not give 
rise to dynamical instability; they only generate gaps in 
the spectra. Based on this general conclusion, we have 
analyzed dynamical instability for two extreme limits, 
»<1 and c<t). 

The limit, v -C 1, corresponds to the cases shown in 
the first column of Fig. [3] For this limit, we can approx- 
imate the phonon spectrum and the anti-phonon spec- 
trum with the ones given in Eq. I|5.12l) . By equating them, 
£+(q— 1) = £_ (q), for the degeneracy, we find that the dy- 
namical instability should occur on the following curves 

k = ^q 2 c + q^/4 + y/(q - l) 2 c + (q - l) 4 /4 . (5.13) 

These curves are plotted as solid dots in Fig.0 and they 
fall right in the middle of the thin dark-shaded areas. To 
some extent, one can regard these thin dark-shaded areas 
as broadening of the curves Ij5.13|l ■ It is noted in Ref . |2l| 
that the relation (|5.13|) is also the result of e+(q — 1) + 
£+(— q) = 0, which involves only the physical phonons. 
Therefore, the physical meaning of Eq. (|5.13|l is that one 
can excite a pair of phonons with total energy zero and 
with total momentum equal to a reciprocal wave number 
of the lattice. 

The other limit, c <C v, is the cases plotted the first 
row of Fig.|SJ The open circles along the left edges of 
these dark-shaded areas are given by 

Ex(k + q)- Ex(k) = Bi(fc) - E 1 {k - q) (5.14) 



where E\{k) is the lowest Bloch band of 
1 d 2 

Ho = --^j + vcos(x). (5.15) 

In this linear periodic system, the excitation spectrum 
(phonon or anti-phonon) just corresponds to transitions 
from the Bloch states of energy E\(k) to other Bloch 
states of energy E n (k + q), or vise versa. The above 
equation is just the resonance condition between such 
excitations in the lowest band (n — 1). Alternatively, we 
can write the resonance condition as 

EAkj+E^k) = E ± (k + q) + E ± (k - q) . (5.16) 

So, this condition may be viewed as the energy and mo- 
mentum conservation for two particles interacting and 
decaying into two different Bloch states E\{k + q) and 
E\(k — q). This is the same physical picture behind 
Ea. (j5~T3f . 

One common feature of all the diagrams in Fig. [5] is 
that there is a critical Bloch wave number kd beyond 
which the Bloch waves <f>k are dynamically unstable. The 
onset instability at kd always corresponds to q — 1/2. In 
other words, if we drive the Bloch state 4>k from k = 
to k = 1/2 the first unstable mode appearing is always 
q = ±1/2, which represents period doubling. Only for 
k > kd can longer wavelength instabilities occur. The 
growth of these unstable modes drives the system far 
away from the Bloch state and spontaneously breaks the 
translational symmetry of the system. 

The dynamical instability of the lower loop states is 
studied in Ref. |2l| . The range of k, where dynamical 
instability occurs, is found decreasing with increasing c 
and decreasing v. 

D. Tight-binding approximation 

We have so far presented a theory of how to determine 
the Landau and dynamical instabilities for the system of 
BECs in optical lattices. It is valid for the full ranges of 
parameters c and v as long as the mean-field treatment 
is adequate. It is, nevertheless, worthwhile to consider 
the system in certain extreme limits, where analytical 
results are available. The particular limit that we now 
look into is the tight- binding limit, where the lattice is so 
strong that the BEC system can be considered as a chain 
of trapped BECs that are weakly linked. This limit was 
studied in Ref.[24| with its focus on dynamical instability. 
In this case, the BEC trapped in each well of the lattice 
can be approximately described by 

tp n (x) = tp°(x-2mr), (5.17) 

where the real wave function ip {x) is very localized in 
the well [0,27r] such that 

i- J dx<p n (x)<p m (x)*5 mn . (5.18) 
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The tight-binding approximation is to write 

^(x,t)=^2W n (t)(p n (x). (5.19) 

n 

Substituting it into the GP equatio n ( 1 2.3(1 . we find after 
neglecting a trivial energy constant |24j. 



.dW n 
''~dt~ 



-t(W n 



+i 



U\W n \ 2 W n , (5.20) 



where 



, . 1 dlfn dip n+1 
2 dx dx 



+ tp n tp n+ ivcos(x)) , (5.21) 



and 



U = 



2tt 



dxifi{x) 



(5.22) 



No surprise, the tight-binding equation l|5.20|l has also a 
Bloch wave solution 



W% — exp (i2kmr - fit) 



(5.23) 



with fi(k) = U — 2icos(2fc7r). The instabilities of this 
Bloch wave solution can be determined by following pre- 
cisely the method described in the previous subsections. 
We write down the perturbed state 

W n {t) = W*{t) + {ue 2lqn7T + v*e- 2lqn *y 2kn7 <- tlt . (5.24) 

which leads to an M matrix (similar to Eq.JOJ 



M k (o)-( L{k + q) U 
M *W-y U L{-k + q) 

In the matrix, 

L(k + q) = U + 4tsin(g7r) sin[(q + 2k)ir] 



(5.25) 



(5.26) 



By straightforward calculations, we find the eigenvalues 
Of M k (q) 

^k(q) = + 4isin 2 (g7r)cos(2£;7r)± 



±y/U 2 + [2<sin(2g7r)sin(2fc7r)] 2 , (5.27) 
and the Bogoliubov spectrum (eigenvalues of a z Mk(q)) 
Ek(q) = 2t sin(2q7r) sin(2fc7r) ± 2sin(g7r) x 

2tC/cos(2fc7r) + At 2 sin 2 (g7r) cos 2 (2/ctt) . (5.28) 

Based on these two results, one can easily determine the 
stability phase diagram of the Bloch wave solution l|5.23|l . 
which is plotted in FigEl When U/2t < 1, there is a 
right border line for dynamical instability (see the top 
right corner of Fig|Hfa)); it is given by 



sin 2 (g7r) 



U 



2£|cos(2fc7r)| ' 



*4 



(5.29) 



This border disappears for U/2t > 1. For any value of 
U/2t, the left border of dynamical instability is at k = 
1 /4 for any q while the boundary of Landau instability 
is described by 



cos 2 (q7r) = cos(2fc7r)^^ + cos(2fc7r)^ 



k < 



1 



(5.30) 



Note that the border line k = 1/4 is consistent with 
our analysis in the last subsection. It is the result of 
Eg. (|5. 1411 when one uses the tight- binding band energy 
E(k) oc cos(2fc7r). 

It is evident that the phase diagram FigEl is a nat- 
ural extrapolation of the results shown in Fig[Sl and it 
also follows the trend in Fig[S]as one increases the inter- 
action: the overall area of instability decreases while the 
dynamical instability spreads more in the area of Landau 
instability. 



cr 1/4 



cr 1/4 




FIG. 6: Stability phase diagram in the tight-binding limit. 
Two different typical cases, U/2t < 1 and U/2t > 1, are 
shown in (a) and (b), respectively. 



E. Experiments 

This dynamical instability was observed in a recent 
experiment 12], whose setup is schematically drawn in 
Figt3 A cigar-shaped BEC was formed in a magnetic 
trap superimposed with an optical lattice. This essen- 
tially prepared a BEC in a Bloch state at k = 0, the 
center of the Brillouin zone. The magnetic trap was kept 
on to keep the quasi-one dimensional configuration, but 
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its center was suddenly shifted by Ax along the longitu- 
dinal direction. This is equivalent to displacing the whole 
BEC off the center of the harmonic trap then releasing 
it. 



magnetic trap 




optical lattice 



FIG. 7: Schematic drawing of the experimental setup in 
Refs.EHII. 



According to the semi-classical theory [31 



dE(k) 
dk 



-iv 2 x . 



(5.31) 



where E{k) is the system energy and lo is the longitudi- 
nal frequency of the magnetic trap. The BEC begins to 
oscillate in both the real space and the momentum space 
(Brillouin zone). Since E{k) oc k 2 near the center of the 
Brillouin zone, both oscillations are harmonic when the 
displacement Aa; is small, as observed in the experiment. 
However, when Aa: was over a certain critical value, say, 
Aa; c , it was observed that the oscillations were disrupted 
and became dissipative. 

This dissipative behavior is caused mainly by the dy- 
namical instability in our opinion although the authors 
of Ref . 01 attributed it to Landau instability. When the 
displacement Ax is small, the oscillation is confined near 
the center of the Brillouin zone, where Bloch states are 
stable as shown in Fig.[S] With the increase of Aa;, the 
Bloch state will be driven into the dangerous zone of dy- 
namical instability; the oscillation then can be disrupted 
by the growth of modes of instability. In our units, this 
experiment has v ~ 0.2 and c ~ 0.02, which is in the 
regime where the dynamical instability is rampant. Fur- 
ther experiments in Ref.|l2j found that Aa; c increases 
with the decrease of the lattice strength v, and there was 
no dissipation when the density of BEC was low. These 
two observations basically rule out Landau instability as 
the main cause of dissipation. As we see in each row 
of Fig.0 the Landau instability, as represented by the 
light-shaded areas, is affected very little by the change 
of v, thus Aa; c should see no changes within experimen- 
tal error when v changes. Also clear in Fig.|Sl when c 
is small, almost all Bloch states have Landau instabil- 
ity. As a result, more severe dissipation would have been 
observed for smaller c if the Landau instability was the 
main contributor. 

Very recently, there was another experiment with 
a similar setup, but the optical lattice used is much 



stronger [23. The advantage of strong optical lattices is 
that dynamical instability is more severe, rendering Lan- 
dau instability less relevant. As a result, the interpreta- 
tion of the experimental observation is more definite and 
much less controversial than the one in Ref. 0] . 

There are other possible ways to observe the dynami- 
cal instability. One way is to prepare a Bloch state then 
monitor how its periodicity changes over with Br agg scat- 
tering of a probing laser beam by the BEC cloud |47j ■ An- 
other way is indirectly through the breakdown of Bloch 
oscillations [ill [ill l32l| . Within proper regime of param- 
eters, the dynamical instability can be very severe and 
thus drive the system far away from a Bloch state within 
a Bloch period. This leads to the breakdown of Bloch os- 
cillations; the observation of this breakdown is certainly 
an observation of dynamical instability. However, one 
must be careful since the nonlinear Landau-Zener tun- 
neling due to the loop structure, as discussed in the last 
section, can also leads to the destruction of Bloch oscil- 
lations. This ambiguity can be avoided by using dilute 
BEC such that c < v, where the loop does not exist. 



F. Discussion on dynamical instability 

As we emphasized above, dynamical instability hap- 
pens only when both the optical lattice and the interac- 
tion between atoms exist. However, it does not imply 
that it is special to BECs in optical lattices. To the 
contrary, it is very general. For example, it has been 
identified in a BEC confined in a torus. In the following, 
we will demonstrate this point with a one-mode problem. 
To provide a slightly different angle, we will present it in 
the language of operator algebra and use the Bogoliubov 
approach^- A discussion on a similar two-mode problem 
can be found in Ref. . 

The one-mode problem is described by the Hamilto- 
nian written as 



H = Aa^a + — a T a T + 



B 



B 



(5.32) 



where A and B are two positive real numbers. The two 
operator and d are generation and annihilation op- 
erators, respectively, satisfying the boson commutation 
relations 



[a, a 



f ] = l, [&,&] = [tf,tf]=0 



(5.33) 



This kind of quadratic Hamiltonians are often found 
as the result of application of perturbation theory to 
a degenerate Bose gas near a condensed state, such as 
Ea. (|5.3() . For this simple one- mode case, the M matrix 
similar to Ea. (|5.4|) is 



M = 



A B 
B A 



(5.34) 



The Bogoliubov approach is to diagonalize the Hamil- 
tonian l|5.32|l with the introduction of a new set of oper- 
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ators, 

a = ub + v*P , a"<=u*tf + vb (5.35) 
For b and to be boson operators, it requires that 



M 2 -M 2 = l. (5.36) 



If u and v satisfy 



o,M\ u \=e( u 
vl V v 



(5.37) 



the Hamiltonian l|5.32jl becomes diagonal in b and w , up 
to a constant, 



(5.38) 



The phonon energy e is easily found by solving 
Ea. H5.37j l. and it is 



enters a new unknown phase and demands a new mathe- 
matical description other than the mean-field theory, the 
GP equation. In Ref . plj . the new phase was called "in- 
sulator" . There may be some experimental indications to 
call it "insulator" ; theoretically, we are yet to develop a 
theory accurately depicting the new phase. Our theory, 
along with the one in Ref. |24L |25| , can only predict when 
or where the dynamical instability sets in and the system 
enters a new phase; it does not give any indication how 
the new phase should look like and behave. 

The dynamical instability is closely related to Lan- 
dau instability. As already briefly stated above, Bloch 
waves with dynamical instability must also have Landau 
instability; on the other hand, Bloch waves with Lan- 
dau instability are not necessarily dynamically unstable. 
This comes from a rigorous result on the relation between 
eigenvalues of M^{q) and <r z Mk{q) as proved in the Ap- 
pendix: when all the eigenvalues of Mk(q) are all positive, 
the eigenvalues £&(<?) of a z Mk(q) are all real. 



B 2 



(5.39) 



When A < B, this energy becomes imaginary; we, again, 
encounter dynamical instability. The simplicity of this 
one-mode problem allows us look more carefully into this 
instability. We notice from Ea. l|5.38jl that H = when 
e becomes imaginary. At the same time, we find with 



Eci. (j5~5?j> that 



0, violating the condition 



(I5.36J1 for b and w to be boson operators. In other words, 
the Hamiltonian (|5.32|) defies the diagonalization by the 
Bogoliubov approach when A < B. 

The physics lies behind the failure of the Bogoliubov 
diagonalization is the dynamical instability: the system 
is intrinsically unstable and, therefore, does not have any 
phonon excitation. We can look into this in a different 
representation. We write a and a* in terms of space and 
momentum operators, 



a=-j=(x + ip), a f = -^=(x - ip) . 



(5.40) 



As a result, the Hamiltonian assumes the standard form 
of a harmonic oscillator, 



H = 



A + B 



x 2 



A-B 



-P 



(5.41) 



VI. SUMMARY 

In summary, we have examined the superfluidity and 
many of its related properties of BECs in optical lattices. 
We studied the tunneling between Bloch bands and found 
that the tunneling is not zero even in the adiabatic limit 
when the repulsive interaction c is strong enough. This 
dramatic tunneling behavior is due to the loop structure 
in the Bloch band, a result of strong superfluidity We 
investigated Landau instability and dynamical instability 
of the system. With a phase diagram, we showed where 
Landau and dynamical instability occur and how they 
change with the variation of the two system parameters, c 
and v. We have discussed the experimental observation of 
dynamical instability and its implication in this system. 
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It is evident here that, when A < B, the system is a 
harmonic oscillator with a negative mass, an unstable 
system. 

A quantum system with dynamical instability, such as 
a harmonic oscillator with a negative mass, does not exist 
in nature. Its appearance in the theoretical study only 
signals the onset of instability of the underlying unper- 
turbed states. In our case, these unperturbed states are 
the BEC Bloch states. Once these Bloch states develop 
dynamical instability, it means that the GP equation, 
which predicts the existence of these Bloch states, is no 
longer a good description of the system. The system 



APPENDIX A: MATHEMATICAL PROPERTIES 

OF M k {q) AND a z M k {q) 

For simplicity, we will drop subscript k and refer to 
these two matrices simply as M(q) and a z M(q) in this 
appendix. We are interested in the properties of their 
eigenvalues as defined in 



and 



M(q)X(q) = X(q)X(q), 



a z M(q)Y(q) = e(q)Y(q) . 



(Al) 



(A.2) 
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The eigenvalues X(q) are all real since the matrix M(q) is 
Hermitian. On the other hand, the matrix a z M(q) is not 
Hermitian and its eigenvalues e(q) are not necessarily 
all real. However, due to its special structure, the 
eigenvalues e(q) enjoy some very interesting properties. 
There is also a very important relation between the two 
sets eigenvalues, X(q) and e(q). 

(i)If £ is one eigenvalue, then e* is also an eigenvalue 
of o~ z M(q). In other words, if the matrix a z M(q) has 
complex eigenvalues, they appear in conjugate pairs. 

This becomes evident when we write the matrix in a 
form, where all its elements are real. Note that, simi- 
lar to the Hamiltonian operator in the left hand side of 
Ea. (|3.1|) . a z M(q) is periodic in space so that one can ex- 
press it in the space of Fourier series. The expression 
is 



a z M(q) = 




(A.3) 



where 

A mn = ^ 1^ dxe- imx C(k + q)e mx , (A4) 
27r Jo 

B mn = ^- r dxe- imx 4>l{ X y nx , (A5) 
2tt Jo 

and 

C mn = ±- dxe~ mx C(-k + q)e mx . (A6) 
^ Jo 

One can check by straightforward calculations that 
all the elements, A mn , B rnn , and C mn are real. This 



particular expression is very convenient for numerical 
calculation: after neglecting Fourier terms of orders 
higher than a cut-off order, say, N, one reduces a z M(q) 
to a real (AN + 2) x (4iV + 2) matrix. To our experience, 
N = 10 is sufficient. 

(ii) If e is an eigenvalue of a z M(q), then — e* is an 
eigenvalue of o~ z M(—q). 

Write Eq.(A.2) in another form, 

(:)-(::)■ <-> 

By some simple algebraic manipulations, one imme- 
diately observes that — e* is an eigenvalue of a z M(—q) 
with the corresponding eigenvector being {v* , u*}. 

(iii) When all the eigenvalues X(q) are positive, the 
eigenvalues e(q) must be all real.. 

We multiply both sides of Eq.(A.2) by a z then by Y\ 
and obtain 

Y^M(q)Y = e{q)Y^a z Y . (A8) 

Since the matrix M(q) is Hermitian, the left hand side of 
the above equation is always real. As a result, if an eigen- 
value e(q) is complex, then for its corresponding eigen- 
vector Y, one has Y*M(q)Y = Y^a z Y = 0. 

When all eigenvalues of M(q) are positive, the left hand 
side of Eq.(A.8) must be positive. In this particular case, 
one can conclude that the eigenvalues e(q) are all real and 
share the sign of Y'a z Y. 
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